! Two dimensinal Riemann problem solution
! ref. 1.  John D. Anderson. Fundamentals of Aerodynamics, 5th ed.. New York: McGraw-Hill, 2011. Chp. 9
!      2. John D. Anderson. Hypersonic and High-Temperature Gas Dynamics, 2nd ed.. Washington D. C.: AIAA, 2006. Chp. 2
program riemann()
   implicit none
   integer::NPM
   real*8::MaA,pA,thetaA,rhoA,TA,VA
   real*8::MaB,pB,thetaB,rhoB,TB,VB
   real*8::MaC,pC,thetaC,rhoC,TC,VC
   real*8::MaD,pD,thetaD,rhoD,TD,VD
   real*8::pi,gamma,theta,e,ep,thetaE,temp0
   
   NPM=11   !number of Prandtl-Meyer expansion waves
   pi=3.141592653
   gamma=1.4
   e=1e-3
   ep=1e-2
   
   MaA=6.3
   pA=16000
   thetaA=-1*pi/180
   TA=1100
   rhoA=pA/(287.058*TA)
   MaD=6.5
   pD=25000
   thetaD=1*pi/180
   TD=1200
   rhoD=pD/(287.058*TD)
   
   if (pA < pD) then   ! A-->B shock, D-->C expansion, thetaA < thetaD
      theta=thetaD-thetaA+e
      call shock(gamma,theta,MaA,pA,rhoA,TA,VA,MaB,pB,rhoB,TB,VB)
      if (pB > pD) then	!double-shock
         theta=e
         call shock(gamma,theta,MaA,pA,rhoA,TA,VA,MaB,pB,rhoB,TB,VB)
         thetaB=theta+thetaA
         thetaE=abs(thetaB-thetaD)
         call shock(gamma,thetaE,MaD,pD,rhoD,TD,VD,MaC,pC,rhoC,TC,VC)
!         write(*,*) pB,pC,thetaB*180/pi
         do while (abs((pB-pC)/pD) > ep)
            theta=theta+e
            call shock(gamma,theta,MaA,pA,rhoA,TA,VA,MaB,pB,rhoB,TB,VB)
            thetaB=theta+thetaA
            thetaE=abs(thetaB-thetaD)
            call shock(gamma,thetaE,MaD,pD,rhoD,TD,VD,MaC,pC,rhoC,TC,VC)
!            write(*,*) pB,pC,thetaB*180/pi
         end do
         write(*,*) 'double shock'
         write(*,*) 'B',MaB,pB,rhoB,TB,VB,thetaB*180/pi,'deg'
         write(*,*) 'C',MaC,pC,rhoC,TC,VC,thetaB*180/pi,'deg'
      else if (pB == pD) then	!single-shock
         MaC=MaD
         pC=pD
         rhoC=rhoD
         TC=TD
         VC=VD
         thetaB=thetaD
         write(*,*) 'single shock'
         write(*,*) 'B',MaB,pB,rhoB,TB,VB,thetaB*180/pi,'deg'
         write(*,*) 'C',MaC,pC,rhoC,TC,VC,thetaB*180/pi,'deg'
      else   !one shock-one expansion
         thetaB=thetaD+e
         thetaE=abs(thetaB-thetaD)
         call PME(NPM,gamma,thetaE,MaD,pD,rhoD,TD,VD,MaC,pC,rhoC,TC,VC)
!         write(*,*) pB,pC,thetaB*180/pi
         do while(abs((pB-pC)/pD) > ep)
            theta=theta+e
            call shock(gamma,theta,MaA,pA,rhoA,TA,VA,MaB,pB,rhoB,TB,VB)
            thetaB=theta+thetaA
            thetaE=abs(thetaB-thetaD)
            call PME(NPM,gamma,thetaE,MaD,pD,rhoD,TD,VD,MaC,pC,rhoC,TC,VC)
!            write(*,*) pB,pC,thetaB*180/pi
         end do
         write(*,*) 'one shock-one expansion'
         write(*,*) 'B',MaB,pB,rhoB,TB,VB,thetaB*180/pi,'deg'
         write(*,*) 'C',MaC,pC,rhoC,TC,VC,thetaB*180/pi,'deg'
      end if
   else if (pA > pD) then   ! A-->B expansion, D-->C shock, thetaA < thetaD
!      write(*,*) 'Please exchange region A and D, and be sure pA < pD'
      theta=abs(thetaA-thetaD-e)
      call shock(gamma,theta,MaD,pD,rhoD,TD,VD,MaC,pC,rhoC,TC,VC)
      if (pC > pA) then	!double-shock
         theta=e
         call shock(gamma,theta,MaD,pD,rhoD,TD,VD,MaC,pC,rhoC,TC,VC)
         thetaC=-theta+thetaD
         thetaE=abs(thetaC-thetaA)
         call shock(gamma,thetaE,MaA,pA,rhoA,TA,VA,MaB,pB,rhoB,TB,VB)
!         write(*,*) pB,pC,thetaC*180/pi
         do while (abs((pC-pB)/pA) > ep)
            theta=abs(-theta-e)
            call shock(gamma,theta,MaD,pD,rhoD,TD,VD,MaC,pC,rhoC,TC,VC)
            thetaC=-theta+thetaD
            thetaE=abs(thetaC-thetaA)
            call shock(gamma,thetaE,MaA,pA,rhoA,TA,VA,MaB,pB,rhoB,TB,VB)
!            write(*,*) pB,pC,thetaC*180/pi
         end do
         write(*,*) 'double shock'
         write(*,*) 'B',MaB,pB,rhoB,TB,VB,thetaC*180/pi,'deg'
         write(*,*) 'C',MaC,pC,rhoC,TC,VC,thetaC*180/pi,'deg'
      else if (pC == pA) then	!single-shock
         MaB=MaA
         pB=pA
         rhoB=rhoA
         TB=TA
         VB=VA
         thetaC=thetaA
         write(*,*) 'single shock'
         write(*,*) 'B',MaB,pB,rhoB,TB,VB,thetaC*180/pi,'deg'
         write(*,*) 'C',MaC,pC,rhoC,TC,VC,thetaC*180/pi,'deg'
      else   !one shock-one expansion
         thetaC=thetaA-e
         thetaE=abs(thetaC-thetaA)
         call PME(NPM,gamma,thetaE,MaA,pA,rhoA,TA,VA,MaB,pB,rhoB,TB,VB)
!         write(*,*) pB,pC,thetaB*180/pi
         do while(abs((pC-pB)/pA) > ep)
            theta=abs(-theta-e)
            call shock(gamma,theta,MaD,pD,rhoD,TD,VD,MaC,pC,rhoC,TC,VC)
            thetaC=-theta+thetaD
            thetaE=abs(thetaC-thetaA)
            call PME(NPM,gamma,thetaE,MaA,pA,rhoA,TA,VA,MaB,pB,rhoB,TB,VB)
!            write(*,*) pB,pC,thetaC*180/pi
         end do
         write(*,*) 'one shock-one expansion'
         write(*,*) 'B',MaB,pB,rhoB,TB,VB,thetaC*180/pi,'deg'
         write(*,*) 'C',MaC,pC,rhoC,TC,VC,thetaC*180/pi,'deg'
      end if
   else if (thetaA <= thetaD) then  !pA == pD
      VA=MaA*sqrt(gamma*287.058*TA)
      VD=MaD*sqrt(gamma*287.058*TD)
      thetaB=(thetaA+thetaD)/2   !double equal shock
      theta=thetaB-thetaA
      call shock(gamma,theta,MaA,pA,rhoA,TA,VA,MaB,pB,rhoB,TB,VB)
      MaC=MaB
      pC=pB
      rhoC=rhoB
      TC=TB
      VC=VB
      write(*,*) 'double equal shock'
      write(*,*) 'B',MaB,pB,rhoB,TB,VB,thetaB*180/pi,'deg'
      write(*,*) 'C',MaC,pC,rhoC,TC,VC,thetaB*180/pi,'deg'
   else
      write(*,*) 'No Riemann problem encountered!'
   end if
   
end program riemann

!shock
subroutine shock(gamma,theta,Ma1,p1,rho1,T1,V1,Ma2,p2,rho2,T2,V2)
   implicit none
   real*8::gamma,theta,Ma1,p1,rho1,T1,V1,Ma2,p2,rho2,T2,V2
   real*8::beta,Mn1,temp0,temp1
   V1=Ma1*sqrt(gamma*287.058*T1)
   if (theta == 0) then
      Ma2=Ma1
      p2=p1
      rho2=rho2
   	T2=T1
   	V2=V1
   	return
   end if
   beta=theta+1e-3
   Mn1=Ma1*sin(beta)
   temp0=tan(beta-theta)/tan(beta)
   temp1=(gamma+1)*Mn1*Mn1
   temp1=(2+(gamma-1)*Mn1*Mn1)/t1
   do while(abs((temp0-temp1)/temp0) > 1e-2)
      beta=beta+1e-3
      Mn1=Ma1*sin(beta)
      temp0=(temp0-temp1)/temp0
      temp0=tan(beta-theta)/tan(beta)
      temp1=(gamma+1)*Mn1*Mn1
      temp1=(2+(gamma-1)*Mn1*Mn1)/temp1
   end do
!   write(*,*) beta*180/3.141592653
   Ma2=gamma*Mn1*Mn1-(gamma-1)/2
   Ma2=(1+(gamma-1)*Mn1*Mn1/2)/Ma2   !Mn2**2
   Ma2=sqrt(Ma2)/sin(beta-theta)
   p2=p1*(1+2*gamma*(Mn1*Mn1-1)/(gamma+1))
   rho2=rho1*(((gamma+1)*Mn1*Mn1)/(2+(gamma-1)*Mn1*Mn1))
   T2=T1*(p2/p1)/(rho2/rho1)
   V2=Ma2*sqrt(gamma*287.058*T2)
end subroutine shock

!Prandtl-Meyer Expansion
subroutine PME(NPM,gamma,thetaE,Ma1,p1,rho1,T1,V1,Ma2,p2,rho2,T2,V2)
   implicit none
   integer::NPM,i
   real*8::gamma,thetaE,Ma1,p1,rho1,T1,V1,Ma2,p2,rho2,T2,V2
   real*8::thetaEi,niu,niu1,niu2,Pt,Tt,temp
   V1=Ma1*sqrt(gamma*287.058*T1)
   thetaEi=thetaE/NPM
   temp=Ma1
   do i=1,NPM
      Ma2=temp+1e-4
      niu1=niu(gamma,temp)
      niu2=niu(gamma,Ma2)
      do while(abs((niu2-niu1-thetaEi)/thetaE) > 5e-3)
         Ma2=Ma2+1e-4
         niu1=(niu2-niu1-thetaEi)/thetaE
         niu1=niu(gamma,temp)
         niu2=niu(gamma,Ma2)
      end do
      temp=Ma2
!      write(*,*) Ma2,thetaEi*i*180/3.141592653,thetaE*180/3.141592653
   end do
   Pt=p1*(1+(gamma-1)*Ma1*Ma1/2)**(gamma/(gamma-1))
   p2=(1+(gamma-1)*Ma2*Ma2/2)**(gamma/(gamma-1))
   p2=Pt/p2
   Tt=T1*(1+(gamma-1)*Ma1*Ma1/2)
   T2=(1+(gamma-1)*Ma2*Ma2/2)
   T2=Tt/T2
   rho2=p2/(287.058*T2)
   V2=Ma2*sqrt(gamma*287.058*T2)
end subroutine PME

real*8 function niu(gamma,Ma)
   implicit none
   real*8::gamma,Ma
   real*8::temp0,temp1
   
   temp0=(gamma+1)/(gamma-1)
   temp1=Ma*Ma-1
   niu=sqrt(temp0)*atan(sqrt(temp1/temp0))-atan(sqrt(temp1))
end function niu
